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We apply the postquasistatic approximation, an iterative method for the evolution of self- 
gravitating spheres of matter, to study the evolution of dissipative and electrically charged dis- 
tributions in General Relativity. The numerical implementation of our approach leads to a solver 
which is globally second-order convergent. We evolve nonadiabatic distributions assuming an equa- 
tion of state that accounts for the anisotropy induced by the electric charge. Dissipation is described 
by streaming out or diffusion approximations. We match the interior solution, in noncomoving co- 
ordinates, with the Vaidya-Reissner-Nordstrom exterior solution. Two models are considered: i) a 
Schwarzschild-like shell in the diffusion limit; ii) a Schwarzschild-like interior in the free streaming 
limit. These toy models tell us something about the nature of the dissipative and electrically charged 
collapse. Diffusion stabilizes the gravitational collapse producing a spherical shell whose contraction 
is halted in a short characteristic hydrodynamic time. The streaming out radiation provides a more 
efficient mechanism for emission of energy, redistributing the electric charge on the whole sphere, 
while the distribution collapses indefinitely with a longer hydrodynamic time scale. 

PACS numbers: 04.40.Nr,04.25.-b,04.25D- 



cr 



> 

o^ 
O 
(N 

O 
O 

> 

X 



I. INTRODUCTION 



Despite their apparent simplicity, 1-f-l models of the 
fluid dynamics of compact objects in numerical relativity 
can include realistic transport mechanisms and equations 
of state. Renewed interest on electric charge in stars has 
driven the numerical integration of the Einstein-Maxwell 
(EM) system. Current integrators of the EM system are 
in comoving coordinates jT], and seem to be limited to 
one dimensional numerical solvers [2] , [5] of the May and 
White family [4] , [5] . Because of the obvious interest in 
three-dimensional situations, it is desirable to use non- 
comoving coordinates. Numerical simulations to explore 
the relevance of electric charge in the process of dissipa- 
tive and anisotropic (viscous) gravitational collapse are 
desirable as well. 

The numerical solution of Einstein equations in 3-1-1 
dimensions is an essential tool for the investigation of 
strong field scenarios of astrophysical interest (see [6J and 
references therein). Numerical relativity has led to the 
discovery of critical phenomena in gravitational collapse 
[7] , allowed the study of binary black holes and neutron 
stars [8j-[TT] and the development of relativistic hydro- 
dynamics solvers |12j . among other major achievements. 
The main limitation currently faced by realistic models 
in numerical relativity is the computational demand in 
three-dimensional evolution. Computationally less in- 
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tensive ID models still remain an interesting alternative 
and help narrow the search in the parameter space for 
general solvers. These simplified systems provide a nec- 
essary test bed to study the phenomena expected in fully 
realistic three-dimensional configurations. 

In this paper, we study a self-gravitating spherical dis- 
tribution of charged matter containing a dissipative fluid. 
We use noncomoving coordinates and follow the method 
reported in [14] . Herrera et al. realized that this method 
was equivalent to going one step further from the qua- 
sistatic regime, and consequently has been named the 
postquasistatic approximation (PQSA) after [15]. 

The essence of the PQSA was first proposed in [T5] 
using radiative Bondi coordinates and it has been exten- 
sively used by Herrera and collaborators [T7]-[2^. In the 
context of charged distributions of matter the original 
method was used as well ^\, [Mj, [3S]. The approach 
is based on the introduction of a set of conveniently de- 
fined "effective" variables, which are the effective pres- 
sure and energy density, and a heuristic ansatz on the 
latter [H]. By QSA we mean that the effective vari- 
ables coincide with the corresponding physical variables 
(pressure and energy density). In Bondi coordinates the 
notion of QSA is not evident: the system goes directly 
from static to postquasistatic evolution. In an adiabatic 
and slow evolution we can catch-up that phase, clearly 
seen in noncomoving coordinates; this can be achieved 
using Schwarzschild coordinates [TS]. If the configura- 
tion is leaving equilibrium, the PQSA description seems 
to be enough and can be used as a test bed in numerical 
relativity [26]. Its systematic use of local Minkowskian 
and comoving observers, named Bondians, was used to 



reveal a central equation of state in adiabatic scenarios 
[27] . and to couple matter with radiation [28j . 

In self-gravitating systems the electric charge is be- 
lieved to be constrained by the fact that the resulting 
electric field should not exceed the critical field for pair 
creation, 10^^ V cm~^ [53]. This restriction in the critical 
field has been questioned [30]-I3S] and does not apply to 
phases of intense dynamical activity with time scales of 
the order of (or even smaller than) the hydrostatic time 
scale, and for which the QSA [33] is clearly not reliable as 
in the collapse of very massive stars or the quick collapse 
phase preceding neutron star formation (see |35j and ref- 
erences therein). Electric charge has been also studied 
mostly under static conditions [5S], [5S], [37], [3S]. It is 
of recent interest the charged quasiblack holes [39] and 
its extension to quasispherical realization [40]. Distribu- 
tions electrically charged can be considered in practice as 
anisotropic [5T], [32] ■ Authors combine anisotropy and 
electric charge [3S], [35], [33] but not as a single entity 
by means of an equation of state. 

The electric field has been postulated to be very high 
in strange stars with quark matter [3S], [H], although 
other authors suggest that strange stars wouldn't need 
a large electrical field [371. The effects of dissipation, 
in both limiting cases of radiative transport, within the 
context of the QSA, have been studied in [iS]. Using this 
approximation is very sensible because the hydrostatic 
time scale is very small, compared with stellar lifetimes, 
for many phases of the life of a star. It is of the order of 27 
minutes for the sun, 4.5 seconds for a white dwarf and 104 
seconds for a neutron star of one solar mass and 10 Km 
radius [35]-[nD]. However, such an approximation does 
not apply to the very dynamic phases mentioned before. 
In those cases it is mandatory to take into account terms 
which describe departure from equilibrium, i.e. a full 
dynamic description has to be used 51 . 

In this paper we consider that the electric charge can 
be seen as anisotropy |41], but not any anisotropy as we 
shall see. For certain density ranges, locally anisotropic 
pressure can be physically justified in self-gravitating 
systems, since different kinds of physical phenomena may 
take place, giving rise to local anisotropy and in turn re- 
laxing the upper limits imposed on the maximum value 
of the surface gravitational potential [HI] . The influence 
of local anisotropy in general relativity has been studied 
mostly under static conditions (see |53] and references 
therein; [54j ) . Herrera et al. [HS] have reported a general 
study for spherically symmetric dissipative anisotropic 
fluids with emphasis on the relationship among the Weyl 
tensor, the shear tensor, the anisotropy of the pressure 
and the density inhomogeneity. 

On the other hand, it is well known that different 
energy-momentum tensors can lead to the same space- 
time [SB]-[BD]. For instance, viscosity can be considered 
as a special case of anisotropy [61]. Here we illustrate 
this idea for the Einstein-Maxwell system under spher- 
ical symmetry. To accomplish that program we use the 
total energy characterization as in [3S] and [55] . The elec- 



tric energy (or pressure) contributes to the fluid in a such 
way that the electrically charged perfect fluid is equiva- 
lent to an anisotropic fluid under certain conditions. 



Massive stars evolve emitting massless particles (pho- 
tons and/or neutrinos). Neutrino emission seems to be 
the only plausible mechanism to carry away the bulk of 
binding energy of a collapsing star, leading to a black 
hole or neutron star [53]. It seems clear that the free- 
streaming process is associated with the initial stages of 
the collapse, while the diffusion approximation becomes 
valid toward the final stages. Observation from super- 
nova 1987A indicates that the regime of radiation trans- 
port prevailing during the emission process is closer to 
the diffusion approximation than to the free streaming 
[64] . Heat flow is usually considered as proportional to 
the gradient of temperature. This assumption is very 
sensible because the mean free path of particles respon- 
sible for the propagation of energy in stellar interiors is 
very small as compared with the typical length of the 
object [23]. 



Although some transport equations in the relaxation 
time approximation have been proposed (see for instance 
[65] and references therein), the evolution of temperature 
profiles in the context of general relativity remains an un- 
solved problem. Therefore, we avoid stating any explicit 
evolution equation for heat fiow in this investigation. Re- 
cently, some progress have been achieved by Herrera and 
collaborators on the study of dissipation via an appro- 
priate causal procedure (see for example [55] |B6 68 ) . In 
the present investigation we obtain the zeroth order level 
of approximation for heat flow profiles, which will serve 
as the basis of a future investigation which includes dissi- 
pation in a realistic way using the Miiller-Israel-Stewart 
theory [5M7^ . We already see an important advantage 
using the PQSA studying configurations with anisotropy 
(induced by shear viscosity) , streaming out and heat flow 
processes in spherical collapse: an observer using radiat- 
ing coordinates to study heat flow misses some important 
details. 



The paper is organized as follows. In Sec. II we write 
the field equations for Bondian observers to show how 
the electric charge induces anisotropy, matching with the 
exterior Reissner-Nordstrom-Vaidya solution, and write 
the surface equation, following the PQSA protocol [l5] . 
[73]. In Sec. HI we present a summary of the nu- 
merical methods employed. In Sec. IV we show local 
and global tests of numerical convergence and illustrate 
the PQSA integration procedure with two nonadiabatic 
charged models. In Sec. V we conclude with some re- 
marks. 



II. THE EINSTEIN-MAXWELL SYSTEM 

A. Field equations for Bondian observers 

To write the Einstein field equations we use the hue 
element in Schwarzschild~like coordinates 



ds' 



e^dt^ - e^dr^ 



(1) 



where v — i^(t, r) and A — X(t,r), with {t,r,9,(j)) = 
(0,1,2,3). In order to get physical input we introduce 
the Minkowski coordinates (r, x, y, z) by [73] 

dr = e'^/'^dt, dx = e^/'^dr, dy = rdO, dz = rsinM^. (2) 

In these expressions i' and A are constants, because they 
have only local values. Next we assume that, for an ob- 
server moving relative to these coordinates with velocity 
u! in the radial direction, the space contains a nonstatic 
distribution of matter which is spherically symmetric and 
consists of charged fluid of energy density p, pressure p, 
electric energy density pe, radiation energy flux q diffus- 
ing in the radial direction, and unpolarized radiation of 
energy density e. Thus, the energy-momentum tensor is 



T„ 



{p+p)u^u^-pgf_,^+eli_,l^ + q^u^ + q^Uf_, + E^^, (3) 



where Uq,, Zq,, g^ are the 4- velocity, the 4-null vector 
and the heat flux 4-vector respectively, which satisfy 
u"ua = 1, qau" = 0, l"la = 0, and Ef^^, is the elec- 
tromagnetic energy-momentum tensor 



E, 



fii' 



F "F 



o F F'^'^ 



(4) 



is naturally interpreted as the charge within the radius 
r at the time t. We define the function s{t,r) by F*^ = 

gg-(A+,.)/2/^2^ with 



s(t,r)= / 47rrV*e~(^+'')/2^r. 



(9) 



The conservation of charge inside a sphere comoving 
with the fluid is expressed as 



u"s,a = 0. 



(10) 



We can write the conservation equation in a more suitable 
form for numerical purposes 



dr 

S,t + -T:S,r = 0, 

dt 



(11) 



where the velocity of matter in the Schwarzschild coor- 
dinates is 



dr 
di 



= c.e('^-^)/2. 



(12) 



The contravariant components of the 4-velocity, the 
heat flux 4-vector and the null outgoing vector are 



-1^/2 



we 



-A/2 



(l-Cj2)l/2 * (l-Cj2)l/2 






e 



-A/2, 



1 (1-^2)1/2°* + (1_ ^2)1/2"- 



and 



If" = e-^/Hi" + e-^/^St!. 



(13) 

(14) 
(15) 



We write the field equations for the Einstein-Maxwell 
system in relativistic units (G = c = 1) as follows: 



where F^^ is the Maxwell field tensor, which satisfies the 
Maxwell equations: 



and 



-P[/il/;<T] — 



iV^F^n,. = 47^y^J^ 



(5) 



(6) 



where the semicolon (;) and the comma (,) represent co- 
variant derivative and partial differentiation with respect 
to the indicated coordinate, respectively; J'* = uu^ is 
electric current 4-vector and a the electric conductivity. 
Because of the spherical symmetry, only the radial elec- 
tric field i^*'" = —F^*^ is nonzero. On the other hand, the 
inhomogeneous Maxwell equations become 




P + 



87rr4 



1 2 

---{e~^[2i^rr + l^^r ~ ^,r'^,r + -{v,r - A.; 

327r ' r 



(16) 



(17) 



2 7*^5(1'+^) 



and 



= 47rr J e2 



s., = -4^rVe^(-+^), 



(7) 



(8) 



where J* and J^ are the temporal and radial components 
of the current 4-vector, respectively. The function s{t^r) 



and 



-e-''[2X,tt + XA\t-yM 



uj . , (l + a;2) 1 



l-UJ 

87rr 



l-u;2 



l-uo 



(18) 



(19) 



Anisotropy induced by electric charge 



where 



To write the field equations in a form equivalent to an 
anisotropic fluid, we introduce 



1 - 2^/r, 



(20) 



fi{t,r) ^m{t,r) 



2r' 



(21) 



m being the total mass 



Thus the field equations ( 16 )-( 19 1 read 



P 



47rr^ 



(22) 



1 

Sirr 



2^i 2m 



(23) 



IQ-KT I 



2\ (m,,. - ii/r) 



87r(r - 2^) 



M,tt 



3M^t M,ti',t 



(r - 2a.) 



(24) 



and 



47rr 



(1 -2/i/r)3e" 



where the conservative variables are 

'2 2wg 1 



P = 






1 -w2 



l-w' 



5(Pr+P) + 7" 



and the fiux variable 



P 



Pt - P^ 
1-^2 



2ujq 



(25) 



(26) 



(27) 



.^, (28) 

1 — ui 



as in the standard ADM 3+1 formulation. Within the 
PQSA p and p are referred as effective density and effec- 
tive pressure, respectively. 



Equations (22)-(25l are formally the same as for an 
anisotropic fluid, with p = p + pe, Pr — p — Pe, Pt = 



p + Pe and the electric energy density p^ — E 



where 



E ~ s/r is the local electric field intensity. If we define 
the degree of local anisotropy induced by charge as A = 
Pt—Pr — 2/Oej the electric charge determines such a degree 
at any point. 



From ( 22 ) and ( 25 1 we easily obtain 



dp 
~di 



( fir 

_4^r2 ^ j^ ^ g^ ^ ^)j (^ _ 2/Vr)i/2e''/2 

(_ at 



This equation, known as the momentum constraint in the 
ADM 3+1 formulation, expresses the power across any 
moving spherical shell. 



It can be shown that 

{p + p){4Trr^p + p) 



P.I 



r{r — 2/i) r 



A'!Tr{r — 2p) 



Pm 



ip - Pt) 



P.tv.t 



r-2p 



(30) 



This equation is the same as for an anisotropic fluid [75] 
and is a generalization of the hydrostatic support equa- 
tion, that is, the Tolman-Oppenheimer-Volkoff (TOV) 
equation. Equation ( 30 1 is equivalent to the equation of 



motion for the fluid in conservative form in the standard 
ADM 3+1 formulation [55]. Equation ([30| leads to the 
third equation at the surface (see next section) ; up to this 
point is completely general within spherical symmetry. 

To close this sub-section we have to mention that we 
assume the following equation of state (EoS) [75] for 
nonadiabatic modeling and only as initial-boundary da- 
tum: 



Pt -Pr 
where C is a constant. 



C{p + ij){4Trr^p + p) 
{r ~ 2m) 



(31) 



C. Junction conditions 

We describe the exterior spacetime by the Reissner- 
Nordstrom-Vaidya metric 



ds\ 



1 



2M{u) , e , ^ , 



r r 

r^(de^ + sin^ dd(f>^ 



du + 2du dr 



(32) 



where Ai{u) is the total mass and £ the total charge, 
and u is the retarded time. The exterior and interior 
solutions are separated by the surface r — a{t). In order 
to match both regions on this surface we use the Darmois 
junction conditions. Demanding the continuity of the 
first fundamental form, we obtain 



1 



2M 
a 






and 



main) =M 



= -Xa- 



(33) 



(34) 



(35) 



(36) 



The subscript a indicates that the quantity is evaluated 
at the surface. In this work, we will use the continuity of 
the independent components of the energy-momentum 
flow instead of the second fundamental form, which have 
been shown to be equivalent [34- and it is simpler to apply 
to the present case. This last condition guarantees the 
absence of singular behaviors on the surface. It is easy 
to check that 



Pa =qa, 



(37) 



which expresses the continuity of the radial pressure 
across the boundary of the distribution r = a(t). 



D. Surface equations 
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FIG. 1: Anisotropic parameter h as a, function of the total 
electric charge £, calculated using ( |31[ ) evaluated at the surface 
for the Schwarzschild-like model. The initial conditions are 
a(0) = 5.0, m„(0) = 1.0, a;<,(0) = -10"^ 
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Following the protocol sketched in TF we write the 
surface equations evaluating (12 1, (29) and (30) at the 
surface of the distribution. The first and second surface 
equations read 



da 



d\Ji.a 

dt 



1- 



2Aia 



= -L 



e- da 
2a2 di' 



with 



L=[Q + E{l+OJa)] 1- 



2a*q 



(38) 



(39) 



(40) 



where E = Ana^ea and Q = Aira^qa- 

We need a third surface equation to specify the dy- 
namics completely for any set of initial conditions and a 
given luminosity profile L{t). For this purpose we can use 
the field Eq. (18) or the conservation Eq. (30) written 



in terms of the effective variables, which is clearly model 
dependent. 



FIG. 2: Initial profile of the charge function s, calculated 



using (31 1 for the Schwarzschild-like model. The initial con- 
ditions are a(0) = 5.0, ma{Q) = 1.0, u]a{Q) = -10"^, with 
£ = 0.5, which corresponds to h = 0.8966 



III. NUMERICAL METHODS 

Once the surface equations are integrated using a stan- 
dard method, as Runge-Kutta of 4th order (RK4), we 
have to integrate the conservation equation (fTTl) to ob- 
tain all the physical variables inside the source. Thus the 
conservation equation 



■s,t 



dr 
di' 



(41) 



is a wave-like equation that can be integrated nu- 
merically using the Lax method (with the appropriate 
Courant-Friedrichs-Levy (CFL) condition). The evolu- 
tion of the conservation equation is restricted by the sur- 
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FIG. 3: Initial profile of the energy density p (a varia- 
tion relative to the outermost value, multiplied by 10^) for 
the Schwarzschild-like model I. The initial conditions are 
a(0) = 5.0, niaiO) = 1.0, uja{0) = -10"^ with I = 0.5, which 
corresponds to h = 0.8966. 



FIG. 5: Initial profile of the matter velocity dr/dt (multiplied 
by 10"^) for the Schwarzschild-like model I. The initial con- 
ditions are a(0) = 5.0, ma{0) = 1.0, uja{0) = -10"^, with 
^ = 0.5, which corresponds to h — 0.8966. 
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FIG. 4: Initial profile of the radial pressure p (multiplied by 
10^) for the Schwarzschild-like model I. The initial conditions 
are a(0) = 5.0, ma(0) = 1.0, ujaiO) = -10"^ with £ = 0.5, 
which corresponds to h = 0.8966. 



FIG. 6: Initial profile of the heat flow q (multiplied by 10®) 
for the Schwarzschild-like model I. The initial conditions are 



a(0) = 5.0, ma{0) = 1.0, a;„(0) = 
corresponds to ^ = 0.8966. 



-10" 



with £ — 0.5, which 



face evolution and is implemented as follows 



o"-|-l 



1 

2'^ 



j+i 



^j^i) 



5^ 
25r 



dry 



=j+i 



^3-1 J 



(42) 
The superscript n indicates the hypersurface t — nSt and 
the subscript j the spatial position for a comoving ob- 
server at r = j6r. In order to integrate the conservation 
equation we need to specify a boundary-initial condition. 
The PQSA is a seminumerical method where the radial 
dependence is determined from a static interior solution 
and kept the same during the evolution. The problem 
is typically reduced to a system of ordinary differential 
equations (ODEs) at the surface of the distribution of 
matter. This system is integrated in time using the RK4 
method. Therefore we can calculate exactly any physi- 



cal variable at the interior. In our specific case, for the 
Einstein-Maxwell system, the approach requires addi- 
tionally the integration of the conservation of the elec- 
tric charge (Eq. (41)) using the Lax method (Eq. (42)). 
The conservation equation is an evolution equation con- 
strained by the system of ODEs at the surface. Thus, 
the numerical convergence of the whole algorithm must 
be of 2nd order accuracy. 

The implemented algorithm at the surface (basically 
that of RK4) for the specific models was verified with 
satisfaction from a physical point of view and within a 
reasonable numerical error (see subsection IV. A). If the 
numerical solution for the electric charge function is sta- 
ble and globally convergent to 2nd order, as shown in 
subsection IV. B, the problem surely is well-posed [77j . 
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FIG. 7: Evolution of radius a for the Schwarzschild-like model 
I. The initial conditions are a(0) = 5.0, ma(0) = 1.0, u!a{0) = 



-10" 



with h = 0.90 and /i = 0.96, corresponding to £ = 



0.491 and £ = 0.314, respectively. 



FIG. 9: Evolution of the matter velocity dr/dt at the surface 
for the Schwarzschild-like model I. The initial conditions are 



a(0) = 5.0, ma{0) = 1.0, w„(0) = 
corresponds to ft = 0.8966 



-10" 



with £ — 0.5, which 



1.955 




& 1.925 

i-i 

(a 

a 

" 1.915 



1.905 



2.0 





5 


00 




4 


95 




4 


90 


■H 
■H 
Xi 


4 
4 
4 


85 
80 
75 




4 


70 




4 


65 




4 


60 



, , . 










\ 


\ 


\ 



12 3 



4 5 
Time 



9 10 



FIG. 8: Evolution of the energy density p (multiplied by 10^) 
at the surface for the Schwarzschild-like model I. The initial 
conditions are a(0) = 5.0, ma{0) = 1.0, uja{0) = — lO"'^, with 
£ — 0.5, which corresponds to h — 0.8966 



FIG. 10: Evolution of radius a for the Schwarzschild-like 
model II. The initial conditions are a(0) = 5.0, ma{0) = 1.0, 
'^a(O) = — 10~^, with h = 0.98 and h = 1 corresponding to 
£ = 0.222 and £ = 0, respectively. 



IV. TESTING AND MODELING 

To illustrate the method let us consider a 
Schwarzschild-like model. Following the protocol 
for the PQSA [TS], [73] the interior solution has the 
effective density 



fit), 



(43) 



where / is an arbitrary function of time. The expression 
for p is 



P+i/5 



1 



-pr 



h/2 



k{t), 



(44) 



P + P V 3 
where fc is a function of t to be defined from the boundary 



condition ( 37 ) , which now reads, in terms of the effective 



variables, as 



Pa = pa^l + (qa + ea.)(l + ^a)^ - (1 + ^a) 



87ra4' 



Thus, (44) and (45) give 






3\ V5C-Xs(l-2Ma/a)''/2 J' 



with 



e = 



1 [r/a) 

a 



h/2 



(45) 



(46) 



(47) 
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FIG. 11: Evolution of the energy density p (multiplied by 10^) 
for the Schwarzschild-like model II. The initial conditions are 



a(0) = 5.0, ma(0) = 1.0, uja{0) 
corresponds to h = 0.9839. 
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FIG. 13: Evolution of the matter velocity dr/dt for the 
Schwarzschild-like model II. The initial conditions are a(0) = 



5.0, m„(0) = 1.0, LOa{Q) 
sponds to ft = 0.9839. 
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FIG. 12: Evolution of the radial pressure p (multiplied by 10^) 
for the Schwarzschild-like model II. The initial conditions are 



a(0) = 5.0, ma{0) = 1.0, uja{0) 
corresponds to h = 0.9839. 
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FIG. 14: Evolution of the radiation flux e (multiplied by 10*) 
for the Schwarzschild-like model II. The initial conditions are 
a(0) = 5.0, m„(0) = 1.0, LJa{0) = -10"^, with £ = 0.2, which 
corresponds to h = 0.9839. 



where h — I — 2C 



,2 I -|^ A^g 
a 



Xs = 6iujl + l)'^+2{Q + E){l+uj,f-{l+u:l)-, (48) 






Using ( 22 ) and ( 23 1 it is easy to obtain expressions for fi 
and v: 



/^ = tJ-air/aY 



(50) 



and 



i^s = 2(3w^ + l)^+2(Q+£;)(l+a;,)2-(l+c.2)_. (49) 
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FIG. 15: Evolution of the charge function s for the 
Schwarzschild-like model II. The initial conditions are a(0) = 
5.0, ma{0) = 1.0, uja{0) = -10~^, with £ = 0.2, which corre- 
sponds to h = 0.9839. 
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to construct the initial profile for the charge function. 
From (31) evaluated at the surface we obtain h{i) (see 
Figure ijT Figure 2 displays the charge function s for the 
following set of initial conditions 

a(0) = 5, TOa(0) = l, c^a(O) = -10-3, (53) 

with i = 0.5, which corresponds to h = 0.8966. 

A. Model I 

In the diffusion limit (e = 0) we choose to = 1.0 and 
S = 0.01, Mr = 10-2. We tested that the algorithm 
is correct at the surface (that is, locally) by verifying 
that the pressure is equal to the given gaussian pulse 
with an accuracy of about IQ-^^. This allow us at least 
to be confident about the implemented algorithm at the 
surface. However, as a double-check, from an strictly 
numerical point of view, we show in Table I a proper 
convergence test to 4th order of the RK algorithm at the 
surface. For this test we use the gravitational potential 
at the surface as the required norm, Af — Ha/O'- Thus, it 
can be shown that the rate of convergence is 



n = log2 






(54) 



where J\fc, Nm and A// are values of N for a coarse, 
medium, and fine time step Ai, respectively (scaling as 
4:2:1, [78], [79]). This corresponds to a local convergence 
test for the RK4, giving a convergence rate n w 4, as 
expected. 



t 


A/'c-A/'^(10-^°) 


Nm-Nf[lQ-^^) 


n 


1.0 


-0.022 


-0.142 


3.964 


1.5 


-0.573 


-3.586 


3.998 


2.0 


-1.148 


-7.180 


3.999 



TABLE I: Proper convergence of the surface gravitational po- 
tential; the expected value for n is 4. 



FIG. 16: Initial and final charge function s profile for the 
Schwarzschild-like model II. The initial conditions are a(0) = 
5.0, m„(0) = 1.0, LOa{Q) = -10"^, with i = 0.2, which corre- 
sponds toh = 0.9839. 



which correspond to the Hamiltonian constraint and to 
the polar slicing condition in the ADM 3-1-1 formulation, 
respectively. It is necessary to specify one function of t 
and the initial data. 

We choose L{t) to be a Gaussian 



L = Lne-(*-*«)'/^' 



(52) 



with Lq = Mr/VSTT, which corresponds to a pulse radi- 
ating away a fraction of the initial mass M^. It is easy 



For the initial setting corresponding to Figure 2 the 
interior docs not evolve because the velocity becomes a 
complex value. However, we can display for analysis the 
initial set of physical variables in Figures 3~6. Under 
these conditions the surface evolves without anomalies as 
shown in Figures 7-9. We have looked for an electrically 
charged initial configuration for which all regions evolve; 
we found one with a total electric charge « 10^^. 



B. Model II 

Let us consider now a Schwarzschild-like model in the 
streaming-out limit ((7 = 0). We choose Iq = 5.0 and E = 
0.25, Mr = 10""'^. We tested again that the algorithm is 
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correct at the surface by verifying that the pressure is 
equal to zero with an accuracy of about 10~^^. For this 
model we show the global rate of proper convergence in 
Table II. For that purpose we construct the following 
norm with the electric charge function s 



M^ / s^dr. 
Jo 



(55) 



For instance, we choose grids of 10, 20 and 40 nodes, 



t{w-') 


.^c(10-3) 


M„ (10-^) 


A/"/ (10~') 


n 





3.5423 


3.5322 


3.5297 


2.0152 


1 


5.1726 


5.1564 


5.1524 


2.0398 



TABLE II: Proper convergence of the norm ( |55[ ); the expected 
value for n is 2. 

and a RK4 time step of 10~^ in a proportion of 4:2:1, re- 
spectively. To reach the same monitoring time the CFL 
time step have to be 6t = K5r, where i^ is a constant of 
order of one, in proportion 4:2:1 as well. Consequently, 
the number of RK4 time steps required to apply the for- 
mula (54) is in proportion 1:4:16. The global convergence 
test (RK4 -I- Lax) gives a rate of n w 2, as expected. 

The initial setting shown in Figure 2 is valid for this 
model too, but with £ — 0.2 (corresponding to ft, = 
0.8966). A larger total electric charge does not allow 
an interior evolution because unphysical values develop. 
Figures 10-16 show these results. The collapse is un- 
avoidable after emitting an energy equivalent to 10 % of 
the initial mass, which decreases the energy density (the 
contrary occurs in the diffusion limit |80j). The electric 
charge is redistributed in the whole body. Note that the 
electric charge on each comoving shell starts moving until 
it reaches a stationary state while the whole distribution 
collapses. Observe how the gradient of charge decreases 
and becomes linear with the advance of time. 



V. CONCLUDING REMARKS 

We consider the evolution of a self-gravitating spher- 
ical distribution of charged matter containing a dissi- 
pative fluid. The use of the PQSA with noncomoving 
coordinates allows us to study electrically charged fluid 
spheres in the diffusion and the streaming out limits as 
they just depart from equilibrium. From this point of 
view, the PQSA can also be seen as a nonlinear pertur- 
bative method to test the stability of solutions in equi- 
librium. We have shown that our seminumerical imple- 
mentation is globally second-order convergent. 

Our results indicate that the dissipative transport 
mechanisms and the equation of state chosen to treat 
electric charge as anisotropy are crucial for the out- 
come of gravitational collapse. We want to stress: 
i) the straightforward manner in which we connected 
anisotropy with electric charge using an EoS; ii) how the 
EoS is used in practice as an initial-boundary condition; 



iii) that departing from the same static solutions we find 
very different evolutions. Increasing the amount of total 
charge £ results in lower values of the anisotropy param- 
eter h, approaching zero. The zero limit is not reached 
because the system breaks down for some limit value of 
the total electric charge. When the transport mechanism 
is diffusive, it was not possible to establish why the sys- 
tem can be set initially but not evolved, except at the 
surface. From the results, the field equations seem to 
be imposing restrictions within the context of diffusion 
and electric charge (or anisotropy) and permit only bub- 
bles of charged matter. Otherwise the electric charge (or 
anisotropy) has to be very small. In the streaming out 
limit the situation is quite different, as is expected. The 
interior is evolved for a total charge that is 200,000 times 
that of the maximum charge permitted in the diffusion 
limit. Coupling of matter with radiation is not strong 
enough to prevent the collapse and the system efhciently 
radiates a large quantity of mass. The system clearly 
departs equilibrium and collapses. Electric charge con- 
tributes to the collapse in the same way that anisotropy 
with tangential pressure greater than radial pressure fa- 
vors the collapse |73| , irrespective of the transport mecha- 
nism. In any case electric charge has to be huge to change 
the fate of the gravitational collapse. The electric charge 
is redistributed in a such way that its gradient decreases 
toward the surface and becomes unexpectedly linear and 
stationary, with the advance of time. There is a critical 
total electric charge (or anisotropy parameter) for which 
the system evolves constrained by the Einstein-Maxwell 
system of field equations. 

Beyond the models we want to stress some features 
about our framework. First, the luminosity profiles are 
given as Gaussian but they can be provided from obser- 
vational data. To keep physical variables on appropriate 
values in the diffusion approximation, the pulse has to 
be narrow in comparison, while the streaming out limit 
allows for a wider pulse. Second, the EoS used in this 
work is not essential. Other EoS as initial-boundary 
data should fit well, understanding that it represents 
anisotropic matter. Third, from the observational point 
of view, temperature profiles are desirable as input data, 
but they are not available in the PQSA method when 
electric charge is taked in account. 

We considered the dissipation by viscosity and heat 
flow separately [73], [5D], in order to isolate similar ef- 
fects with different mechanisms. In this work we consid- 
ered heat flow/streaming out and anisotropy induced by 
electric charge, pointing to the most realistic numerical 
modeling in this area [35,. The results constitute a defi- 
nite first cut to more general situations using the PQSA, 
including dissipation, anisotropy, electric charge, heat 
fiow, viscosity, radiation fiux, superficial tension, tem- 
perature profiles and study their influence on the gravi- 
tational collapse. Numerical issues apart, the inclusion of 
superficial tension [SlJ together with a highly compressed 
Fermi gas [TS] , and more realistic thermal processes [5S] , 
is of current interest in astrophysics [3F . Cooling times 
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of smooth or crusty surfaces may be the way to differen- 
tiate strange stars from neutron stars [47]. The PQSA 
can be used to model these situations. This investiga- 
tion is an essential part of a long-term project which 
tries to incorporate the Miiller-Israel-Stewart theory for 
dissipation and deviations from spherical symmetry, spe- 
cially when considering electrically charged distributions. 
Besides being interesting in their own right, we believe 
that spherically symmetric fluid models are useful as a 
test bed for more general solvers in numerical relativity 



[571 [55] • A general three-dimensional code must also be 
able to reproduce situations closer to equilibrium. 
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